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, Sparse grids are tailored to the approximation of smooth high-dimensional functions. 

^ ' On a d-dimensional cube, the number of grid points is TV = 0{h~-^ \ \ogh\'^^^) with a 

mesh size parameter h. The so-called combination technique, based on hierarchical de- 
^Sj ' composition, facilitates the numerical solution of partial differential equations on these 

grids. Key to the convergence analysis are specific multivariate error expansions, which 
we derive in a generic way for linear difference schemes through an error correction 
. , technique employing semi-discretisations. We obtain error formulae of the structure 

^ ■ e — ©(/i^l log/i|'^~^) and illustrate the convergence by numerical examples. 
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1.1 Background 



High-dimensional models play an important role in various applications. A prominent example 
from finance are option pricing problems, which typically involve computation of the expectation 
of the pay-off with respect to a number of underlying stochastic factors (for instance equities or 
interest rates). Commonly used ways to accomplish this are numerical integration through deter- 
, ministic cubature rules or stochastic simulation, or the solution of the corresponding (potentially 

$H ' high-dimensional, if the number of components is large) Feynman-Kac PDE. 

In quantum-mechanics, high-dimensional eigenvalue problems for the Schrodinger equation gov- 
ern the density functions of the states of large molecular systems. More examples with high- 
dimensional state spaces are found in population dynamics or data mining, to name but a few. 

Classical grid based methods on a d-dimensional cube, with mesh width h, say, suffer from the 
curse of dimensionality: the number of unknowns N required to achieve a prescribed error e by 
an order p method grows exponentially like 

N{e) = 0{s-'^/P). 

The analysis in [2] shows that for a hierarchical Tensor product basis the contribution of single 
piecewise polynomial basis functions (to the approximation of sufficiently smooth functions) is 
proportional to the measure of their support. Optimal approximation for a given number of 
degrees of freedom is attained for so-called sparse grids. Their relevance for the solution of PDEs 
was first revealed by Zenger in [14] , and subsequently error bounds for finite clement methods for 
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elliptic problems were derived in detail in [TJ [2] . Recently, optimal convergence rates were shown 
for parabolic equations under very weak assumptions on the regularity of the initial condition for 
a sparse wavelet method and hp discontinuous Galerkin time stepping in |13j . 

In contrast to Galerkin-type methods, the combination technique - first introduced in [7] - 
decomposes the solution into contributions from Tensor product grids. This has the practical 
advantage that only numerical approximations on relatively small conventional grids need to be 
computed. These can be obtained independently and are superposed subsequently, which lends it- 
self to very efficient parallel implementations. This concept has been successfully used in a number 
of applications, e. g. computational fluid dynamics 8 , quantum mechanics [1] and computational 
finance [12]. 

Theoretical results for this extrapolation scheme, however, which inevitably rely on the expan- 
sion of the hierarchical surplus in terms of the grid sizes, have so far only been obtained for simple 
model problems. Bungartz et al. [S] obtain error bounds in terms of the Fourier coefficients for a 
central difference scheme for the two-dimensional Laplace equation. Pflaum f^|TD] employs Sobolev 
space techniques for the combination solution of a finite element method and prove asymptotic er- 
rors of the form log/i|~^ in the L2-norm for general linear elliptic equations in two dimensions, 
and first order convergence in the H^-norm. A similar result is also shown for the Poisson equation 
in higher dimensions. Common to these two approaches is the use of semi-discrete solutions to 
derive expansions of the hierarchical surplus. In the first paper Fourier representations for semi- 
discrete solutions to the Laplace problem are introduced as stepping stone from the continuous to 
the discrete version. This corresponds to variational formulations with semi-discrete conforming 
subspaces in [21 [TD] ■ 

We introduce a new framework to derive error bounds for general difference schemes in arbitrary 
dimensions. As in the aforementioned articles, again semi-discrete problems will play an important 
role, but in the somewhat different setting of an error correction scheme that determines a suitable 
error expansion by a simple generic recursion. This will be detailed with the central difference 
stencil for the Poisson equation and with upwinding for the advection equation. For these two 
examples sparse grid error bounds are then derived in terms of the dimension of the problem and 
the smoothness of the solution (measured in C'^-Norms). It will become evident how this extends 
to more general linear elliptic and parabolic equations. We conclude with a thorough analysis of 
carefully chosen numerical examples. In particular, we highlight the dependence of the error on 
the smoothness of the solution and the problem dimension. 

1.2 Hierarchical surplus and combination technique 

Consider the c?-dimensional unit cube = [0, 1]'' and a Cartesian grid with mesh sizes hi = 2^'' 
(corresponding to a level U £ Nq) in direction i — 1, . . . , d. 

For a vector h ~ (/ij)i<i<d denote by Uh the approximation of a function on this grid with 
points x/i = {ijhj) o < < Nj.i < j < d i ~ V^i ^ . Hereby u/i is the discrete vector of 
(approximated) values at the grid points, which is extended to by a suitable interpolation 
operator X as Uh Xu/j. We will ultimately be concerned with a setting where Uh is the finite 
difference solution to a scalar PDF, in which case Uh is different from the vector obtained by 
evaluating the exact solution to the PDF, u, at the grid points. Therefore we denote the latter by 

We now define the family U of solutions corresponding to these grids (see Fig. [T]) by J7 = 
iUii))ieNf, with 

U{i) :=M2-., 

i. e. as the family of numerical approximations Uh on tensor product grids with hk = 2^**= . Then 
the hierarchical surplus is the sequence 

6U -.^Si-.-SdU (1) 
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Figure 1: Full grid AI4, sparse grid A^4, and a possible realisation of a dimension adaptive sparse 
grid M ^. 



with 

U{i) ^fc = 0. 

and Gfe the fc-th unit vector. Obviously the difference operators commute. 

Given a sequence of index sets Ain C Nq, the approximation on level n is now defined as 



If this series converges in absolute terms, the error 

hn~u\\^\\ E su{i)\\< E ll-^^WII 

in a suitable norm || • || is therefore bounded by the remaining surpluses. The choice of an optimal 
refinement strategy is thus thrown back to the control of the surplus. If this is done a posteriori, 
dimension adaptive schemes are obtained. This requires that hierarchical surpluses on finer levels 
are estimated from coarser levels. A heuristic example for integration problems is found in [S]. 
Conversely, an a priori analysis requires analytic estimates for the surplus in terms of the grid 
sizes hi, . . . ,hd and is the focus of this paper. 

In [7] , a splitting into lower-dimensional contributions of the form 

d 

u-Uh^Y^ E 7ji,...j„(-;/iji,---,/ij™)/i^i (4) 

"1 = 1 ill, ■■-,3m} 
C {1 d) 

was proposed to analyse the combination technique. The '•' stands for the arguments of u. As 
motivation, before we turn to the derivation of such expansions, we briefly sketch the immediate 
consequences for the combination technique. 
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If we view u as a sequence that is constant over all refinement levels, and therefore set Su = 0, 
then under assumption Q 

= (5(C/(i) - = O (2-Pl'l) (5) 

with |i| = |i|i := X]fc=i *fei because all lower order terms cancel out through the difference operator. 
Thus an asymptotically optimal choice for the index set is 

Mn - {i e < : |i| < n}. 

The number of elements of Mn is then given through 

the number of nodes in the grid corresponding to an index i is bounded by J^^^j^ (2''= + 1) = O (2'''). 
For this choice of grid we conclude from ([5]) 

\\u - uj < J2 ^ E Oi2-P\'\) - 0(n^-i2-f"). (6) 

ii^Mn |i|>n 

The significance of the result ^ is that although the number of degrees of freedom is only 

Ndof = O (n'^-i2") 

compared to 2'*" on the full grid, the error only deteriorates by a factor of order n'^^^ compared 
to the full grid result 2"^". 

We illustrate the combination formula by a two-dimensional example: 

uq — u(0, 0) 

ui = [it(l,0)-u(0,0)] + [u(0,l)-it(0,0)]+uo 

= [m(1,0) + u(0,1)] -m(0,0) 

U2 = [u{2, 0) - 0)] + 1) - ^(l, 0) - m(0, 1) + u(0, 0)] + [u(0, 2) - m(0, 1)] + ui 

= [u{2, 0) + 1) + u(0, 2)] - 0) + u{0, 1)]. 

The general structure in two dimensions is 

n 71 — 1 

Un=Y^ U{1, n - - E - 1 - 0- 

In one dimension obviously the sparse grid is identical to the full grid. In general, for d> 1 and 
n > d — 1, the combined solution can be written as 

S{n) = ^C/(i), 

|i| — n 

where the difference operator Ji, applied d—1 times to the index n, is the one-dimensional version 
of III), JSJ. Evaluating the coefficients explicitly gives 

n 

Un= ^ fln-; S{n) (7) 

l=n-d+l 

with 

a,^,= (^^l)d-i-^( 0<^<d-l. (8) 
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1.3 Notation and outline 



We see that the error analysis fahs into two main parts: Firstly, an expansion like ^ has to 
be shown for the problem at hand. We consider linear PDEs on ~ [0, l]'' and denote by diU 
the partial derivative of u with respect to Xi. For a multi-index a = (ai, . . . ,ad) let D"u = 
. . . d'^'^u. The appropriate function spaces are those with bounded mixed derivatives 

xi := {u: [0,1]'^^M; i:>'^weCo(/'^)V/3<a} 

Xi(K) := {u^Xi: \\dPu\\^ < K^jS < a] 

xi ;= {u: [0,1]'*^M: L>"ue Co(/'^)VQ;e {0,...,A:}''} 

Xt{K) := {ueX^: \\D'^u\\^<Kyae{0,...,kY}. 

By Ca{I'^) we denote the space of continuous functions which arc zero at the boundary, so we 
require for Xi{K) that all mixed derivatives up to order k vanish at the boundary. 

A detailed error representation will be derived for the finite difference discretisation of the 
Poisson problem and the advection equation. For the Poisson problem 

d 

Au = J2dfu^f in/^ (9) 

u = on dl'^, (10) 

we assume sufhciently smooth / such that u E Xf for the solution of ([5]). For the advection 
equation 

d-l 

Ut + Yd,u = Va; e /''"S i e [0,1] (11) 

1=1 

u(a;,0) =uo(x) Va: € /'^"^ (12) 

u{x,t) = ui{x,t) Va; e 3i : a;, = 0, t e [0, 1] (13) 

we take for simplicity of notation the velocity constant and equal to one in each direction, but it 
will be obvious how the result generalises to the case with variable velocity. We can also identify 
t with Xd in a combined space-time formulation to write X^iLi ^i'^ = with boundary condition 
u{x) = g{x) Va; s.t. 3i : Xi = 0. Here uq, ui (or g, respectively) are required to fulfils some 
compatibility conditions such that u € A^ for the solution u of ([TT|) . 
With one-sided differences 

f^k ^k 

(cfc the unit vector in direction k) we consider central differences 5'^hk^khu Poisson equa- 

tion and upwinding with left-sided differences 5'^ and implicit time-stepping for the advection 
equation. The solution is then extended to I'^ by multi-linear interpolation. For the latter insta- 
tionary problem we will illustrate the difference between a sparse grid in the 'space' coordinates 
only and a 'space-time sparse grid'. 

We write the discretised systems in matrix notation 

A,,u,, = fft, (14) 

with the discrete solution vector u/i, and write u(x^) for the (vector representing the) exact 
solution at the grid points x^, such that the truncation error can be written as A^u(x/i) — th- 
in order to split up the truncation error, we define according semi-discrete operators Ajj'^' "'*™' 
in directions ii, . . . , by 

At-'-'u= E 9fu. (15) 
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for the Poisson problem and likewise for the advection equation. In contrast to the fully discrete 
solutions u;i satisfying , this leads to a system of PDEs 

— 

(with respective boundary conditions), where xi^^^''"'^"^^ are defined on hyper-planes 
I^^i'-*'") ={xel'^ : X,, € {jh,^ : 0<j< hrj^}, I < k < m} 

and f^'i- '*") = ^j^g restriction of / to Denote by Rh := r\^'-'^'> the 

restriction to the grid. Let C'(I'-'^' ' '*'"^) denote the space of continuous functions on these planes 
with maximum norm || ■ ||co and derivatives are defined in directions along the planes. 

The rest of the paper is organised as follows. Section [2] contains the main step of the analysis, 
in which we prove the expansion ^ for the error u(x/i) — u/j at the grid points x/j. We then 
extend the solution to by multilinear interpolation (section [3]) and show that the expansion 
([1]) is maintained pointwise for the error between the exact solution u{x) and the interpolated 
approximation (Xuh){x). The latter result is interesting in its own right as it yields expansions 
for the interpolation error of sufficiently smooth functions on sparse grids and cubature formulae 
that are based on such interpolation properties will be obtained en passent. 

In the second step of the analysis (section H]), these error representations on Cartesian grids are 
used to estimate the error of the combined solution by combinatorial formulae. Exact asymptotic 
expansions in terms of the smoothness of the solution will be given. In section \E\ numerical 
results will illustrate the theory. Section [6] discusses these results and points out future directions, 
extensions to other problem classes and possible applications. 

The author wishes to thank Mike Giles for very illuminating discussions on the subject, in 
particular for pointing him towards the concept of adjoint error correction [H], which proved vital 
for the theory of section [H 

2 Error expansion for finite difference schemes 

The goal of this section is to establish error representations of the form 

d 

u{xh)-Uh^Yl Wju-,j^i^h;hj^,...,hj,JhP^- ...-hP^ (16) 

m=l {Jl , • . • , 3m } 
C {1, - . . , d} 

between the finite difference solution and the exact solution u evaluated on a grid x/j. 
We have in mind linear (possibly degenerate) elliptic equations of the type 

d d 

Au = ttijdidjU + bjdjU + csu = f 

i,j=i j = i 

with X)f j=i "^Mnjaij > 0, and view parabolic equations as a special case, in which we identify the 
'time' coordinate with Xd and set bd — —1, hi — Q for i — — 1, and ai^d = a-d.i = for 

« = 1, . . . , c?. Upwinding in direction d is then equivalent to fully implicit time-stepping. 

It seems instructive to first look at the main principles of the proof, which are generic and 
can be applied to a large class of problems. We restrict this sketch to two space dimensions in 
order to avoid cumbersome notation. We then extend this to higher dimensions and specify the 
requirements for particular examples (Poisson, advection equation), derive the coefficients in the 
expansion (fTH]) and give sharp bounds. 
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2.1 Outline of proof in two dimensions 

Starting point is a consistency assumption of order p of the form 

A,,u ~f^h\- T^^\■,h^) + hi ■ rf (., + Khl ■ Tf^l{-M,h2). (17) 

This win typically be straightforward to see by Taylor expansion. In so doing we are assuming 
that the solution is sufficiently smooth, and we are making such an assumption throughout the 
following considerations. Note that often, e. g. for the Poisson problem, t^^j (•, /ii, ^,2) — 0, but 
this term will be present for mixed derivatives. 

To derive the respective convergence order, one would usually proceed to write 

w(x^) - = hi ■ A-Vf ^(x;,, /ii) + /^^ A-Vf ^(x;,, h^) + h\hl ■ A-\^^}{^h, h,, /la) (18) 

and deduce from the boundedness of A,~^ convergence of the scheme. In this context, however, 
since A^"'^ depends on all grid sizes hi, ... , hd, it is not possible to derive an error expansion of the 
form pop in this way, unless A^^ is known explicitly in terms of the grid sizes and an expansion 
of A.y^^ can be derived. This effectively happens in j3j, where Fourier series of continuous and 
discrete solutions to the Laplace equation are used. Note that this is generally rather cumbersome 
and only possible in special cases where such representations are known. 

This is where the concept of error correction comes into play, and we determine the error terms 
by solving the auxiliary semi-discrete problems 

A««;i(.,/ii) = r(°)(.;/ii) (19) 
A.^u'^w^i-M) = rf(.;/i2) (20) 

on the hyper-planes ij^^^ and respectively (for the definition of these operators see the preceding 
section). Note that wi(-, hi) and W2(', ^2) indeed only depend on hi and ft,2 respectively, because 
they are the solution of a semi-discretisation as in (fTS|) and not of the fully discrete equation. 
Under suitable regularity assumptions, the first terms on the right-hand side of 

Ah{u~ h\wi{-]hi) - hlw2{-;h2)) - f = 

hi (A^^Ki{--hi) - A„u.i(-; hi)) + hi (a.^^'^W2{-, h2) - AhW2{-, /is)) + hlhlr^fi-; hi,h2) 

can be absorbed in the higher order terms by further expanding 

Ai''^ -Ah\wi{-,hi) = hlai.^2{-;hi;h2) (21) 
W2i-;h2) = h\a2;i{-,h2;hi), (22) 



so if we define 
we get 



''1^2 '^i°2 +o-i;2(-;^i;^2) + 0'2;l(-;^2;/jl), 



A,, {u ~ hlwi{-]hi) - hlw2{-;h2)) - / = h^hlT^^^i-; ft.i,/i2). 
Now we can do the final step as we would have done in (fT5|) to obtain 

u(x,i) - Uft = hlwi{xh;hi) + hlw2ixh;h2) + ft-f ft.^wi,2(x;,; /ii, ^.2) (23) 

with 

wi,2(x/i; /ii, ft.2) = A^Vf 2'(x;,;/li,/l2). 
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In higher dimensions, this procedure will be applied inductively. The start of the induction, 
m = 1, is always the consistency assumption (jl7p . The final step to = d + f is always essentially 
equivalent to ([25]) . 

The central proofs in this article go along these lines and it will become clear how this framework 
can be applied to other settings. The main gap that has to be filled concerns the existence of an 
expansion (|2ip . (P^ . and the boundedness of its coefficients, which relies on the smoothness of the 
solutions of p^ . (PO)) . In the following two sections we detail the regularity requirements and give 
explicit bounds on the coefficients in the expansion for the Poisson and the advection equation. 

2.2 Detailed analysis for the Poisson equation 

As an example, we consider central differences for the Poisson problem In this case p — 2. 
2.2.1 Regularity, stability 

Standard elliptic regularity results do not apply here due to possible singularities in the corners. 
Also, we need regularity of the solutions to semi-discrete problems, which are systems of elliptic 
equations of this type, and it is not obvious how estimates independent of the grid sizes can be 
obtained, which is crucial in our context. Therefore we formulate the following lemma. 

Lemma 2.1. Let u be a solution of Au — f with homogeneous Dirichlet data and - j/jg 
solution a|j'^''"'*'"''u|j*^''"'*"'' — f^(*i'---''">) ^yjiili iliQ definitions from section \l.3\) . Then 



Proof. 1. Let |/| < / and v := -\fxi{l-xi). Then A{u + v) = f + f > and hence u + v <0 
(maximum at the boundary), i. e. it < —v < //8. Similarly u > — //8. 

2. Again a semi-discrete maximum principle holds and the analogous result follows by consid- 
ering = J (central differences are exact for quadratic functions). 



Similar estimates for the derivatives of the solution can be obtained by differentiating the equa- 
tion, and using the fact that we are considering function spaces for which derivatives of sufficiently 
high order vanish at the boundaries. 

2.2.2 Consistency 

For completeness we state the expansions (|17l) and ([?!]) ■ ([2^ in detail. Note that the truncation 
error is defined continuously and not just at the grid points. 

Lemma 2.2 (Truncation error of difference stencil). 1. Let u e Xf{K), then 



1- hl|oo<il!/|| 



oo 




lloo<|||/| 



oo 



□ 



d 



k=l 



for some Tk with 



(24) 



forae{QAY, afc = 0. 
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2. Let u e Xi{K) with /3j = 4 V i ^ {ii, . . . , im}, then 



k^{ii,...,i,„} 

again with i24^ , hut now for a G {0,4}'', = Vi G U {fc} and < Pi 

otherwise. 



Proof. Standard Taylor expansion in one variable. 



□ 



2.2.3 Error expansion on grid 

We first prove an error expansion for the finite difference solution at the grid points. 

Theorem 2.3. Let u G Xf{K) he the solution of the Poisson equation and the finite difference 
approximation on a grid x^, . Then 



C {!,... ,d} 



(25) 



lkji,...j™(-;/iji.--->j™)lloo 

Proof. We prove by induction for \ < m < d 
( 



< K- 



rn — 1 



k=l 



5Z ■■■■h'i^'^^^,...,l^{y^h]h^^,...,h^^) (26) 

{il, • ■ ■ , im} 

where Wi^^,,,^i^ (•; /i^j^ , . . . , /i^j,) and r^j....^^^ (•; /i^j , . . . , /ii„) are functions defined on the hyper-planes 
for which the estimates 

P°n,,...,^„(•;/in,•••,/^.,JI|oo < m!8-(™-i)l2-'"if, (27) 
P°^«^l,....^.(•;/l.l,•••,/^^J||oo < fc! S''^ 12-^= if forl<A;<TO-l (28) 

hold if a e {0,4}'' with ai- = 0, 1 < i < (that is along the planes where partial derivatives are 
defined) . 

The case m = 1, 



A,,u(x,,) - fh = YhlTk{xh,hk), 



k=l 



follows from LemmaHIll 1. with \\D°'Tk\\oo < " ^ {0,4}'' with ak = 0. Now assume 

for m > 1 with the bound ([?7|l . Then the solution Wi^^,,,^i^ of 



(ii,...,i,„) 



(29) 



9 



with homogeneous Dirichlet boundary conditions satisfies (Lemnia l2.11 2. and the comments there- 
after) 

P°u;n,...,.™(-;/in,--->.™)lloo <m!8-"12-"if 

with a e {0, 4}'' and — 0, I < j < k and therefore, for m < d, from Lemma [2?2| 2. there exist 
tTii,...,i„;fe(-; /iji, ■ • ■ , hi^;hk) such that 

A|j*'''"''"''wii,...,i^ - A/iWij,..._i^ = ^ ft,fecrii,...,i„;fe(-;^ji,...,/ii„;/ife) 

and 

P°^n,...,^™;fel!oo < m!8-'"12-("+l)if 
for =0, 1 < j <m and = 0. For ni = d obviously Aj^*^' "''™'' — A/, = 0. So define 

ji. • ■ ■ s.t. 

■ • ■ Jm.} U {j} = {il,. . . ,im.+ l} 

and since the sum has m + 1 terms 

for a S {0,4}'' with 0;^^. =0, 1 < j < m + 1. By induction follows for all m and (|25p is 
obtained directly with ni = d. □ 

From the proof of Theorem 12.31 we see immediately the following result: 

Corollary 2.4. The weights ... j-,^ (x^; /ij^ , . . . , /ij^) in i25\) are the restriction of functions 
Wji,...,jm{'] f^jii • ■ • I f^jm) defined on hyper-planes ^ Jqj- -which the bounds 

\\D"w,,_,^{-;h,„...,h,J\\oo < km^'^K 
hold for a e {0,4}'', = Vi G {ii, . . . ,ik}- 

2.3 Advection equation and extensions 

We now turn to the upwind discretisation of equations (|lip - (|13p with fully implicit time-stepping, 
because unconditional stability in the maximum norm is easy to establish here. Analogous results 
to Lemmas 12.11 and 12.21 can be derived and give the following expansion. 

Theorem 2.5. Letu G X2{K) he the solution of the advection equation andxih the finite difference 
approximation. Then 

d 

u(x/,)-u,i=^ ^ Wji,...j„(x,,;/iji,...,/ij„J/ij, (30) 

C {!,..., d} 

with Wj^^,,,j^ e C ^ll^^'i'---^''")^ and 

lkii,...jn.(-;/iii,--->i„)lloo < K—. 

Proof. First derive an expression for the truncation error (the equivalent of Lcmma l2.2p . 

d 

{A - Ah)u = ^ hkTk{-; hk) 

k=l 
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with 

for a e {0, 2}'*, ak = 0, and a stability result (the equivalent of Lemma [2. ip 

< max(||uo||oo, ||ui||oo)- 

The rest follows by the same steps as for the Poisson problem. □ 

We can combine the above results for the Poisson problem and the advection equation to derive 
error formulae and estimates for an upwind discretisation of advection-diffusion equations of the 
form 

b • Vu — Au — cu. 

The truncation error is additive, regularity and stability apply similarly under the assumptions 
made above. The derivation can therefore follow the same steps. 
Instationary problems 

Mt + b • Vu = Am — cu 

or 

d d 

Ut + bidiU — Qidfu — cu 

i=l i=l 

can be seen as a degenerate case with vanishing diffusion in one direction, in which case upwinding 
in this coordinate is equivalent to fully implicit time-stepping (see comment at beginning of the 
section). We can therefore construct a 'space-time' sparse grid, which splits up both space and 
time in a hierarchical way. For other time discretisations, e. g. the second order Crank-Nicolson 
scheme, which has a time step constraint for stability in the maximum norm, time has to be 
treated separately: a sparse grid is used for the space-like coordinates, while the time step has to 
be chosen to satisfy the appropriate stability criterion on the highest space refinement level. 

For non-diagonal diffusion-tensors, it seems difhcult to construct finite difference schemes, for 
which discrete maximum principles hold on anisotropic grids. In this case the above analysis is 
not applicable. We discuss this in further detail in section E) 

More generally, all linear problems that admit smooth solutions and satisfy stability properties 
that can be carried over to the discrete case, fit into this framework. 

3 Multilinear interpolation 

For the definition of an approximation on a sparse grid it is necessary to extend the finite difference 
solution by interpolation. Wc will first show that the interpolation error for a sufficiently smooth 
function by piecewise multilinear splines on a Cartesian grid has an error expansion of the form 
(|3ip . Subsequently we show that the difference between the numerical approximation of the PDF 
(the interpolated finite difference result) and the exact solution still has such an expansion. 

As a central ingredient we first derive a particular partial Taylor expansion. This resembles the 
semi-discretisations of the previous section. 

Lemma 3.1 (Fxpansion of functions with bounded mixed derivatives). Let u G X2, x = {xi, . . . , Xd), 
and introduce a;'^*^ G K'' with jj^.'i = Xj if j S . . - ik}, otherwise, and likewise for 
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s, then 

d d 

u[x) = u(0) + '^^Xidiu[x — x^^'') + ^ XiXj didju{x — x'^^'^^) + 



i # 3 



+ xi- ...■Xddi...ddu{Q) + [ '(^.-s»)5fu(s«)ds. 



=1 







{xi - Si) ■ . . . ■ {xd - Sd) 9i • ■ • 9<iu(s) dsd . . . ds 



Proof. See Appendix [X] □ 

The obvious, but crucial point is that aU terms that are linear in one or more directions are 
interpolated exactly (in these directions). 

Theorem 3.2 (Interpolation of functions with bounded mixed derivatives). Assume u G X| and 
let Xu(Tih) be the multilinear interpolating function on a gridxh- Then 

d 

u{x) ~ {Iu{xh)){x) = Yl Yl aj,^,„j,^{x;hj,,...,hj^)h^.-...-h'^. (31) 



/ , / , "Jl.--- JiTi V-^J "-Jl 1 ■ ■ • J "Jmy-ji ■■■ "■]„ 

m=l {ji,...,jm} 

C {1 d} 



with 
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^ji,...j™(-;^ii,---,/jj„.)l|oo < ( ^ ) \\dj^...dl u\ 



Proof. See Appendix]^ □ 

Remark (Cubature). From approximation results of this form error expansions for cubature for- 
mulae are obtained directly. Since the trapezoidal rule is exact on piecewise multilinear functions, 
the integration error is just the integral of the error terms over [0, l]'^. 

To see that the error expansion is preserved under interpolation of the finite difference 
solution from the grid to [0, 1]'', we split 

u{x) - {IvLh){x) = u{x) - {Iu{y.h)){x) + {Iu{yih)){,x) ~ {I\ih){x) 

= u{x)^{Iu{i^h)){x) + [I{u{xh)-nh)){x). (32) 

The first term u{x) — {Tu{ii.h)){x) is the interpolation error of the exact solution and is given by 
([3l]l . The second term, the interpolation of the discretisation error 

d 

u{-Xh) - Mh = ^ •...•/i2^u;j,,...j,(x,,;/ij,,...,/ijj (33) 

k—l {Jl jm} 

C {1, d} 

on the grid x?i, is problem dependent (see Theorem 12. 3p and needs to be evaluated separately. It 
is straightforward to see that the linear interpolant of the error is bounded by the error at the 
grid points, but it is essential, and more involved, to derive the exact form of the expansion. 
We show this again for the Poisson problem first. 

Theorem 3.3. Let u G Xf{K) be the solution to the Poisson problem and u/,, the numerical 
solution with central differences, Tu^ its multilinear interpolation. Then 

d 

u-IxXh^Y. ^Ji,-,J,^i-^h,„...,hjJhl- ...-h^^^ (34) 

m=l {ji, . . . , Jm } 
C {1, d} 
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with 

lkii,...,i™(-;^ji,---,/ji„)lloo <c K- — , 

and C < 150188 depends on neither the dimension nor the data. 

Proof. See Appendix]^ □ 

Similarly one gets for the advection equation the following result. 

Theorem 3.4. Let u e (i^) be the solution of the advection equation and the numerical 
solution with upwinding and implicit Euler time- stepping. Then 

d 

u-lMh^^ ^ Wj,,...j„(-;/iji,...,/ij„J/iji (35) 

m=l iii, . . . ,jm} 
c {1, d} 

with 

lkji,...,i™(-;^ii,---,/ij,Jlloo <c -K ■ — 

and C < 3/2 depends on neither the dimension nor the data. 

Proof. See Appendix]^ □ 

4 Combined error bounds and asymptotics 

We study now in detail the extrapolation effect that the combination formula has on error terms 
by means of combinatorial relations. This leads to error bounds for the combination solution of 
the order seen in Starting point is the pointwise expansion of the error on tensor product 
grids with mesh sizes h — {hi, . . . , hd), 

d 

™=1 {jl,--.,i,r.} 

(e. g. as shown in the preceding sections for the Poisson and advection equation), where 

K,...,jJ<K yi<m<dy{n,...,j^}ci{l,...,d}. (37) 
Obviously this second step of the analysis is independent of the problem, given that holds. 

4.1 Error bounds 

In [7j Griebel et al. derive from (|36p and ([57]) in the two- and three-dimensional case for p = 2 the 
bounds 

\u -Un\< K2-^" ^1 + (38) 

and 

|.-«„|<if2-2"(^l + |n+|n2), (39) 

respectively. This section is devoted to a generalisation of and ((5^ to error bounds of the 
form 

\u-un\< Kc{d,p)n'^-^2-P'\ 
for arbitrary dimension and arbitrary order. In addition the asymptotic limit 

\U-Un\ 

hm —r-r: 
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will be given. We will use the representation of the combined solution 



w„ = [St'S]{n), (40) 
S{n) = ^{/(i) (41) 

i|— n 

(see section ll.2|) in terms of an iterated application of the one-dimensional difference operator 6i . 
The proof requires a few combinatorial equalities. The following formula for iterated differences 
of products — a discrete version of the product rule for differentiation — will be useful here. 



Proposition 4.1. Let f,g eW^'>. Then 

Si ifg) - E ( • ) sl-'fi'+jKg Vfc e No. (42) 

The straightforward proof (based on induction) is omitted here, but see [TT| . 
Since the number of grids on level n that are involved in the combination solution in dimension 
d is given by 

iV(n,d)=(" + ^^M, (43) 



the following Lemma 14.21 states a necessary condition for the consistency of the combination 
technique (i. e. the sum of coefficients of all grids is 1). 

Lemma 4.2 (Consistency). With N from (J^ Vd e N, Vn e N, n > d - 1 

5f-^N{n,d) = 1. 
Proof. By induction one proves for < k < d — 1 



□ 



The following formula forms a key to the proof of the main result in this section, Theorem 14.41 
To see why, look at n as the sparse grid level and s;2~p' as the error of order p on level I, with 
some coefficient s/ that is collected from the number of grids given by the binomial term. Then 
Lemma [4.31 savs that only the highest order terms 2~p("+'^~i) are left in the sparse grid solution, 
whereas the combination formula cancels out all lower order terms that come from the larger mesh 
sizes on the anisotropic grids. This is a more quantitative version of ([5]). 

Lemma 4.3 (Error representation formula). Let m,d > 1, v : M™ M and for n e No 
F{n):= ^ t;(2-*S...,2-*")2-f*i •...•2-P^™. 

|i| =n 

m-l , s 

sr.^ E «(2-^S---,2-'"). 

|i| = ; 

Proof. See Appendix IB] □ 



Then for d e N 



with 
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After these preparations, the error terms can be estimated conveniently. 

Theorem 4.4 (Error bounds). Assume for all Uh a pointwise error expansion of the form ^36\) 
with \31^ . Then the combination solution \4^U^ fulfils the error estimate 

2K [2^ + 1"-'^^^ 



|u-«n|<^^=^(^^;zrj {n + 2{d~l)f-^2-P-. (44) 
Proof. Because of Lemma 14.21 (the exact solution u is constant over the grid levels) 



and one may write 



The error terms in 



|i|— n ''^ 

u-«„-,5^i ;^(ii-c/(i)). 

|i| — n 



-C^W-E E «,„...,,„(•; 2-'",..., 2-^- )2-P'- .....2- 



C {!,... ,d} 



will now be studied separately with the help of Lemma 14. 3[ applied to 

Fn,...,M ■■= E ^Ji.-,^™(-'2-^- , . . . ,2-^-)2-P'- ..... 2-P'- 

|i| — n 



u-u^=Y. E 5t'Fn,-.M. (45) 

C {!,... ,d} 

The point is that with Lemma 14.31 the factors hi then no longer appear separately, but only in 
their highest order hi ■ . . . ■ hn — 2^". It remains to estimate the coefficients, which leads to the 
polynomial terms in n. From j„ I — follows 

\si\<K[ max\sn+d-i-i\ < K I 

therefore 



m — 1 



1=0 



and with Lemma [4 



^ )(-2)- < Ki )E( "'r' (46) 



n + rf + m — 2 \ \ ^ / TO — 1 



m — 1 



i=0 



\St^Fju-,jM\ < 2-f("+'^-i)(2?' + l)™-iif " + (47) 



(48) 



2P + i\'* 1 if 
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Since the number of terms in is 2^* — 1, we get 



2P + lV ^ 2K 



□ 



Let us study equation for p = 2 and small d. For rf = 1, where the sparse grid is identical 
to the full grid, the estimate reduces to |m — u„| < 2K4~^ (the substitution 2"^ — 1 by 2'* in the 
end of the above proof explains the unnecessary factor 2). In the case d = 2 the leading term is 
given by (|49p . such that in the highest power of n 

5 

\u~Un\ ^ K-n4-", 
in accordance with (|38p . Similarly, in three dimensions one gets 

-1 /5\^ _25 



' ' 2 ^4/ 32 

and recovers p9p . Lower order terms differ due to the numbering of the grids and in fact they 
were not estimated separately from (j48l) onwards. 

More interesting, however, is the dependence of on the dimension d for large d. Keeping n 
fixed, one sees from Stirling's formula, 

kl ~ V2^(- 

that the factor depending explicitly on d grows asymptotically like (5e)'^/Vd. It is an interesting 
and practically very relevant question, whether the bounds in (|44p are sharp asymptotically or the 
exponentially growing constants can be omitted by subtle treatment of the error terms. 

4.2 Improved error bounds and exact asymptotics 

The questions raised at the end of the previous section are addressed in the following corollaries. 
Corollary 4.5 (to Theorem 14. 4p . Under the assumptions of Theorem\4.4\ the sharper bound 



holds for d > 2. 

Proof. See Appendix |B] □ 

Corollary 4.6 (to Theorem 14. 4p . // additionally vi d in I5'6]j is continuous in £ with 

^ '■= i^i,...,d(s 0; • • ■ jO) 7^ 0, then the asymptotic behaviour is 

/gp^lX^^-l d-l 

u-u^ = -v [-^^ (d^2"'" + ^ {n'-'2-n . (50) 
Proof. See Appendix iBl □ 

Remark. In the case v = the leading term proportional to n'^~^2~'''" vanishes. This is for example 
the case if the solution has a low superposition dimension, i. e. is a sum of functions that depend 
only on a subset of the coordinates. 

Remark. ([50)1 shows that for fixed n the coefficients of 2~p" tend to for d cx3. This does 
not mean, however, that the approximation is better in higher dimensions. On the contrary, this 
asymptotic approximation was achieved under the assumption that n is large compared to d and 
the asymptotic range is reached for increasing n in higher dimensions. This behaviour is also 
related to the fact that the number of levels contained in the combination solution grows linearly 
in d. 
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4.3 Sparse grid error bounds for model problems 

Finally we can collect the results from the previous sections to obtain error bounds and asymptotic 
error formulae with explicitly computable constants. 

4.3.1 Poisson equation 

We combine the results from Theorems 14.41 and 13.31 to obtain the following. 



Theorem 4.7. Let u e Xf he a solution of the Poisson problem and Un the sparse grid solution 
on level n. Then 

\u-un\<c- sup WD'^uWoo ■ d ■ •(n + 2(d-l))'^-i-4-" (51) 

|a|^<4 

with c < 121000 and 



u-u,,^5~- n'^-i • 4-" + O {n'i-H-") , (52) 
where c < c||Z?°m||oo with a = (4, . . . , 4). 

The dependence on dimensionality and smoothness becomes apparent. 

4.3.2 Advection equation 

A similar result follows for the advection equation from Theorems 14.41 and 13.41 

Theorem 4.8. Let u £ X2 be a solution of the advection equation and u„ the sparse grid solution 
on level n. Then 

\u-ur,\<c- sup \\D"u\\oo-d- (^) •(n + 2(d-l))'*-i-2-" (53) 

|a|oo<2 V4/ 

with c < 2 and 

M - u„ = 5 • ^ • n'^-^ ■ 2-" + O ^ (54) 

where c < c||_D"u||oo with o; = (2, . . . , 2). 

5 Numerical results 

We illustrate the theoretical findings by numerical experiments, and pay particular attention to 
the asymptotic convergence order, the dependence of the error on the dimension and on the 
smoothness of the solution, as reflected in ([^ and 

5.1 Elliptic problems 

Consider central differences for 

Au = / in [0, 1]'' (55) 
u = g oiid[0,lf. 

We choose the data such that the solution is 

1 JL 

\^{X^-P^f\ (56) 



/ ^ d 

u(x) = exp i-^^K 
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Ai > and p G [0, 1]'', that is f{x) — (~1 + ^iVf) ' u{x). This aUows us to control aU 

derivatives in the hght of (|5ip . In particular 

d 

We choose pi = 0.22081976, = 0.29072005,^3 = 0.28051979, p4 = 0.27032006, ps = 0.24122005, 
Pe = 0.17071947, pr = 0.10101947, ps = 0.09021981 to avoid symmetry effects. 

Another obvious choice for a suitable test case would be a combination of trigonometric func- 
tions. We omit such an example as it basically just reproduces the Fourier analysis of sparse grid 
solutions (see f3|) numerically. 

5.1.1 Convergence order and dimensionality 

We start by looking at the case ([56|) with d > 1 and — 1, i — 1, . . . ,d. Then from above 
max|a(|^<4 ||-D"u||oo — 1 for all d. 

Figure [2] shows a logarithmic plot of the error 

e„ := \u{x^) - Un{x^)\ (57) 

for the sparse grid solution w„ at level n, evaluated at a fixed point x,. Here x^ — (1/2, . . . , 1/2), 
the centre point. 




Figure 2: Pointwise error e„ vs. grid level n (left) and vs. the number of grid points N (right) 
for a sparse grid in dimensions 1 up to 8 (from left to right). The continuous curves are the 
asymptotes as in (1551) , fitted to the data - the estimated parameters are given in Table [TJ 



From the result (|52p we know that 

6„ = 0(n«2-f") 

with q = d — 1 and p = 2. We therefore fit the function 

l{n,p,q,r) ^ —r + q\og2n — pn (58) 

to log2 e„ and determine r, p and q by least squares. The results shown in Table [T] correspond very 
well to the theoretical values. Note that it is very hard to estimate the exponent of the logarithm 
properly. 

We perform the same exercise for the error versus the number of grid points N and fit 

l{N,p,q,r) ^ ~r + qlog^N ~ pN (59) 
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Figure 3: A similar plot to Fig. O but comparing the sparse grid in dimensions 1, 3 and 5 (upper 
graphs in the left plot and lower graphs on the right) to the full grid. The full grid error on a 
given level does not depend significantly on the dimension (three lines on top of each other) , 
whereas the sparse grid requires more refinements to achieve the same accuracy. The right 
plot, however, shows the superior complexity in terms of error reduction per unknown. 



d 


1 


2 


3 


4 


5 


6 


7 


8 


p 


2.000 (2) 


1.938 (2) 


1.905 (2) 


1.970 (2) 


1.871 (2) 


1.901 (2) 


1.944 (2) 


1.982 (2) 


q 


-0.001 (0) 


0.478 (1) 


1.44 (2) 


2.86 (3) 


3.42 (4) 


4.75 (5) 


6.24 (6) 


7.76 (7) 


r 


11.43 


3.97 


3.97 


5.31 


6.04 


8.42 


11.61 


15.28 



Table 1 : Coefficients as in ((55| , fitted to the computed errors by regression. In parenthesis see the 
values from the theory. The corresponding curves are plotted with the data in Fig. [21 



d 


1 


2 


3 


4 


5 


6 


7 


8 


p 


2.000 


-1.889 


1.815 


1.846 


1.715 


1.702 


1.689 


1.668 


q 


-0.002 


2.23 


5.08 


8.74 


10.6 


13.7 


16.7 


19.6 


f 


2.43 


4.75 


9.48 


18.5 


24.9 


35.5 


47.1 


59.2 



Table 2: Coefficients as in ((59|). fitted to the computed errors by regression. 



to the data depicted on the right of Fig.[2]and Fig. [31 The coefficients for the sparse grid are shown 
in Table [H and show that because the asymptotic complexity of the sparse grid is independent of 
the dimension (up to logarithmic factors), the order remains approximately 2. 

For comparison, the lines for the full grid in dimensions 1, 3 and 5 have slope —2 (identically 
to the sparse grid), —0.65 (asymptotically —2/3) and —0.34 (asymptotically —2/5) respectively. 
This is a result of the curse of dimensionality, e — 0{N^'^/'^). 



5.1.2 Dependence on derivatives 

It remains to study the effect of smoothness and (an-)isotropy on the convergence. We consider 
the three-dimensional case and vary Ai, A2 and A3. The leading error term is proportional to 

The evaluation point x* is chosen equal to the point p. Fig. [H shows how the error increases with 
increasing (mixed) derivatives. 

5.2 Parabolic and hyperbolic problems 

We consider 

Ut — i^Au + b • Vu — 
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l0g2 En 



I • 6 8 10 12 14 



■ 


Ai 


= 1, Aa = 1. A;, = 1 


♦ 


Ai 


= 100, A2 = 1. A:, = 1 




Ai 


= 10, Aj = 21), A;, = ,50 




Ai 


= 100, A2 = 100, A:) = 1 




Ai 


= 10000, Aa = 1, A3 = 1 



Figure 4: Pointwise sparse grid error e„ as in ([S7|) on level n for equation with solution (|56p 
for different values for Ai, A2 and A3. 



in two dimensions for different values of v. The case z/ = resembles the advection equation. 
To adjust the smoothness, we choose an initial profile 



u{xi,X2, 0) — < 



(tan^ (f I 



arctan ( tan*^ ( ^ - ) ) r < r < r 



1 r < r 
j r<r 
r > r 



where 



r = (1 - e)r for e g [0, 1] and fc > 1 (see Fig. El left, for e = 0.9). 

The value is 1 in an inner circle with centre (7711,7712) and radius r, and changes to within 
a distance of e. The regions are joined together such that fc — 1 is the order of differentiability 
for e > 0. e = is the discontinuous limit. We choose fc = 5, because from the analysis it 
is known that mixed derivatives of order up to 4 are required for upwinding on a sparse grid. 
Homogeneous Dirichlet conditions u — are set at the boundary. To avoid any effects arising 
from the velocity being aligned with the grid, we choose 61 = 0.31415926535897932385, 62 = 
-0.27182818284590452354, we furthermore set 7711 = 0.5(1 - 61), 7712 = 0.5(1-62) and r = 0.5|b|2, 
such that in the non-diffusive case the profile starts from the upper left quarter of the unit square, 
with the outer circle going through the centre (Fig. [U left, for e = 0.9), and is propagated down 
to the lower right quarter, touching the centre from the other side at T = 1. We evaluate the 
solution at the centre = (1/2, . . . , 1/2) and since the exact solution is unknown, we study the 
surplus 

Cn = |u„+i(a;*) - u„(a;*)|. (60) 

We first consider the case v = 0.1, e = 0.9. The numerical solution is shown in Fig. \5\ First 
order convergence is observed as illustrated by Fig. [B] 



We now take e = and change i' — 0.1, 0.01, 0.001. (Fig. [71 left). Due to the discontinuous initial 
condition, the problem does not fit into the theoretical framework. Nonetheless the smoothing 
property of the diffusion operator suffices to provide sufficient regularity for t > 0. 

Alternatively, taking i' — 0, we let the smooth transition collapse to e = 0.1, 0.01, 0.001 (Fig. [3 
right) . 
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Figure 5: Sparse grid solution at level n = 15 for = 0.1, e = 0.9 at t = (left) and t ~ 1 (right) 
log2 



10 12 14 16 



n 




4 


0.060907 


5 


0.25877 


6 


-4.2474 


7 


0.45337 


8 


2.0783 


9 


-2.8072 


10 


0.76331 


11 


1.2188 


12 


1.4628 


13 


1.5811 


14 


1.6559 


15 


1.7072 


16 


1.7449 


17 


1.7604 


18 


1.8364 



Figure 6: e„ as in ([60|l and convergence rate r„ at refinement level n for ly = 0.1, e = 0.9. 
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= 0.01 
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= 0.001 



logs En 
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10 12 14 16 



* e = 0.1 
■ e = 0.01 
- £ = 0.001 




Figure 7: Convergence of the sparse grid solution for discontinuous initial data for small diffusivity 
(left) and smooth initial data without diffusion (right) for with e„ from ()60p . 



Although the problem is sufficiently smooth, the scales involved are not resolved properly and the 
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Figure 8: Sparse grid solution at level n = 15 for u = 0, e ~ Q a.t t = (left) and t = 1 (right). 



i 



6 8 10 12 14 1.6 18 



n 




4 


2.0381 


5 


1.0974 


6 


0.85633 


7 


-0.61664 


8 


1.7166 


9 


-6.8473 


10 


0.29471 


11 


0.87042 


12 


1.9361 


13 


2.4390 


14 


-5.2997 


15 


-0.64137 


16 


-0.43160 


17 


0.38301 


18 


1.58094 



Figure 9: e„ as in ([HO)) and convergence rate r„ at refinement level n for discontinuous initial data 
without diffusion. 

(mixed) derivatives are too large to reach the asymptotic range at feasible refinement levels. 
Finally, consider the limiting case = 0, e = 0, ie a discontinuous initial profile without diffusion. 



Already in Fig. [8] the problems of the sparse grid become apparent. It is noticeable how the 
anisotropic elements in the sparse grid fail to capture the discontinuous transition at the circum- 
ference of the circle. Fig. [^confirms that convergence (if any) is slow and erratic. 



6 Discussion 

In this article we give explicit error bounds for the sparse grid solution to model problems, and 
at the same time provide a framework that can be applied to more general linear PDEs. The 
ingredients are 

1. sufficiently smooth data to yield solutions with bounded mixed derivatives of required order. 
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2. a discretisation scheme that provides a truncation error of mixed order, 

3. stabiUty of the discretisation scheme, that is a bounded inverse in the maximum norm. 

In most cases the truncation error can be assessed easily by Taylor expansion. The bounded- 
ness of the inverse will be harder to establish, but note that no additional requirements to the 
corresponding full grid case are necessary here. Examples that fall into this category (where dis- 
crete maximum principles are known) are elliptic and parabolic equations. The numerical results 
reproduce the theoretical findings nicely. 

Limitations we encountered concern non-smooth problems, dimensions in excess of eight, and 
non-diagonal diffusion tensors. The latter can often be resolved in practice by diagonalising the 
diffusion tensor and performing a principal component or asymptotic analysis, see e. g. |12j . 
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A Proofs from Section S] 

Proof of Lemma \3.1l Define 

d d 

V := w(a:i, . . . , Xd) — it(0, . . . , 0) — Xi 9iu(a; — x*-*-*) — XiXj didju{x — x'^^'^^) 



- xi ■ . . . ■ Xd di . . . ddu\^^ ^^^0 ^ y2 diu{s^''')dsi - 

- E r ■ ■ ■ r"'"(^n - S.J • ... • (x,,_, - s,,_Jdl . . . ds,,_, ...ds, 

. , Jo Jo 

{21, ... , Id-li 

C{l,....d} 

It is straightforward to see that 

v{xi, . . . ,Xd) = if 3k: Xk^O 

and also 

D°'v{xi, . . . , Xd) = if 3k : Xk — /\ ak = I- 

From 

df...d'dV^o 

one gets v = from which the resuh foUows. □ 

Proof of Theorem \3.SX Without loss of generality consider a point x in the box B [0, hi\x . . .x 
[0, hd\. For all corner points p = (pi, . . . ,pd) E Be := ®i=i{0: ^i} Lemma [XT] to express 

d d 

u{p) = u{x) +^{pi - Xi)diu{x^'-^) + ^ {pi- Xi){pj - Xj)didju{x'^''^^'>) + . . . + 

i=l i.] = l 

i 7^ i 



+ (pi - xi) • . . . • (pd - a;d)5i . . . ddu{x) + E / ~ Sj)9,f u(s^'^) dsj + . . . + 
i 

(pi - si) • . . . • (pd - Sd) 9i . . . dlu(s) dsd... dsi (61) 



Here x = (xi, . . . , x^), and introduce x^^^'---'^''\ gl*!---'*) g with a;^.*i'---'"") — Xj if j G {ii, . . . ik}, 
Pj otherwise, and s^.*^' - ''''^ — gj if j G {ii, . . . ik}, Xj otherwise. 
Then insert in the multilinear approximation 

(iu(x,))(x)= 5: ^^i^.....^^^.«b)=^.(x)+x: E E---.-(P'-) 

C {1 d} 



with 

/ d 



because all the terms in (j6ip that are multilinear are exactly represented. By inserting the points 
p € Be one shows after some calculation the representation 



Pii rp 



yfe + ^fc = 1 

— 1, . . . , m 



(62) 
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with h = {hi-^ , ■ • • , /ii„, ) and 

wy,.{x;h) = / ... / \{{zuh,,-sk) 5?...9^^>--^"))ds™ 

The rest fohows because the maximum of ^^(1 — ^) on [0, 1] is □ 

Proof of Theorem \3.S[ From CoroUarv 12.41 m,-^ iki^h', hi-^ , . . . , hi^) (in ([55]) ) is the restriction of 

a function from hyperplanes - to the grid x^, where for the derivatives in the 

continuous directions 

\\D'-w,,,...^,,{-,h,„...,h,)\\o<^\\DPu\U 

with at = 4, i ^ {ii, . . . , ik} Pi = 4 holds. 

By Theorem 13.21 for points x on the hyper-planes 

{Xwi^^,,,^i^{^h\hi^, . . . ,hi^))[x) = /lii, . . . , /lij + 

lii,...,ik\ji,---,jm (^7 hij^, . . . , hi^ ; /ijj , . . . , hj^)hj_^ ■ . . . ■ h^^ (63) 

{il , ■ • ■ , j,n } 

n{ii, . . . = 

with |7ij_..._i^y^^...j^(x; /iji, . . . , /li^,; /iji, . . . , < ^HL'^wHoo where ai, = 4,q;j„ = 2. Val- 
ues between the hyperplanes are obtained by multilinear interpolation, which introduces no new 
maxima, and thus 



I{u{xh)-Uh) = ^3ll,■■■,^k{■l^^'^,,■■■,h,Jhf^■ ...-hi 



{il,. . . Ak] 
C {!,... 



{il, ■ • ■ 

n{ii, . . . = 

with ||/3ii,...^fcy"i,...,j„ (S ^ii : ■ • ■ ! ^ifc ! ''■Ji j • ■ • 7 ^j,„)l|oo ^ -^96^ 27" 

Putting (f64|) and (I3ip together one gets 



TO \ / 4 \ (to — /)! 



"1=1 {Jl. ■ ■ ■ , Jm} 
C {1, d} 

with 

4™ to! ""^ 

^ (to-Q! m! 1 /384 

< ifJ^e384/27< 150188 i^-^. 
96™ 96™ 

□ 

Proof of Theorem \3.Ji\ Follows by combining Theorems 12.51 and 13.21 similarly to the Poisson case 
by observing that one can write second order terms (from the bilinear interpolation) as first order 
terms by defining 



Pii,...,ik ('! ^ii ! • ■ • I hik)^ii ■ ■ ■ ■ ■ ^ifc ~- Pii,...,ik (' i ^ii i • ■ ■ i hi^^)hi.^ ■ . . . ■ h 
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with |/3ij^....^ij, I < I and the final estimate becomes 

II / ■ , Ml "i! 1 / 8 \ ' Ti! o/,7 3 to! 

w„ , {■■,hi,,...,hi ) oo<K—y 77 — < K—e^'^^ < -K—. 



B Proofs from Section |4] 

Proof of Lemma \4-3\ Inserting gives 



|i|— n 



—pi„ 



= ^ v(2-*\...,2-''")2-P^"=i*'= 
|i|=« 

n 



= E^'2- 



n — l + d — m— 1 
d — TO — 1 



From Lemma FB.ll below. [21 {d ^ d — to, /; — > s/2~p') one obtains 

|i| — n 

For aU j > 

S{2-P'' = {2-P - ly 2-P", 
as one sees inductively (j — > J + 1) by 

^j+i2-pn ^ - ly ^2"P(«+i) _ 2-P") = (2-P - 1)-'"^^ 2" 
Therefore it follows from Proposition 14. II that 



■pn 



m— 1 ^ 

TO — 1 



3=0 
m— 1 



3=0 V / 



2-pin+d-l) (^^ _ 2P + 1)"-! 



i=0 
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Lemma B.l (Differencing formula). Let dGN and f,F G s. t. 

1=0 ^ ' 

1. Then for < k < d 

6'lF{n) = G''{n)+H''{n), 

where 

( k = 



1=0 

2. 



n-l+d-1 
d-k-l 



5iF{n) = fn+d 

Proof. 1. Induction in k: Clearly (fc = 0) 

(5?F(n) = F(n) = H°{n) = G°(n) + H°{n). 
Since F{n) has the form F{n) = X^JLo /(^> '^)) 

n 

5iF{n) = /(n + 1, n + 1) - ^ {f{l, n + 1) - f{l, n)) , 

1=0 

hence 

j-j^, , . , \^ r \ f n-l + d \ f n-l + d-1 \ 

6,F{n) = d-1 )-[ d-1 ) 

n , 
= fn+1 + /i ( 



n-l+d-1 
d 



= G\n)+H\n), 

i. e. the case k — 1. 

For fc > 1 the contributions (A; — » + 1) are 

k 



d-j-1 
k-j 



= /n+fe+l - fn+1 ( _ j + ^ /n+i 

- •'"+^' fc - j + 1 J ■'"+1 fc - 1 



k-j+1 \ k-j 
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and in H^{n) 



— fn+l 



d-1 
k 



+1 



n 

d-l 
k 



n — I + d 
d-k-1 



n-l+d-1 
d-k-2 



n- l+d- 1 
d-k-1 



This completes the result as 

5'l+^F{n) = SiS'^F{n) = Si (G'=(n) + H''{n)) = SiG''{n) + 5iH^{n) 
= G''+\n)+H''+\n). 



2. FromHJ one obtains for A; = d — 1 



1=0 



and 



d-l 



SfF{n) — ^(/n+j+l - fn+]) + fn+l — fn+d- 



Proof of Lemma \4-5\ A binomial of the form (|47p or (j48p , respectively, can be written as 
n + d — 1, fc = m — 1) 



fc 



n(i+7) 



From the inequality between the arithmetic and geometric mean one gets 



furthermore 



1 /■'^ dx 

y - < 1+ / — = 1 + lnfc, 



because the first sum is a lower sum for the integral. This gives 







14 


l + \nk' 




)< 






k \ 



and 



n + d + m — 2 
TO — 1 



< 



1 + ln(TO - 1) 



TO — 1 



m — 1 



The rest follows as in the proof of Theorem 14.41 
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Proof of Lemma \4-.6] For continuous wi....,m we can first show 

lim 



— — = fi,...,m(0, . . . ,0) =: Vq. 

n + m — 1 \ 

m — 1 J 

In the foUowing the index of wi,...,m is omitted for simphcity of notation. Let e > and 



(66) 



Choose no such that 



c := sup \v{hi,...,hm)-vo\. 

0</ll....,/lm<l 



Vzfc >no,l <fc<ni: \v{2-'\ . . . ,2-'"^) - vo\ < - 



and then n sufficiently large such that 
No := 



n — kno + m — 1 
m — 1 



> 1 



71 + m — 1 
m — 1 

V 



The latter is possible, because 



n — kno + m — 1 
TO — 1 



n + TO — 1 
TO — 1 



n 1 



kno 
n + j 



for n — > oo. The idea is now to show that the contribution of the terms that do not lie in an e-ball 
around vq can be neglected. This motivates the splitting 



J2 v{2-^\...,2-^-). 



12k ik = n 
Vfc : 2fc > no 



3k : ik < no 



Because of 



one sees 



E 1= E 



l = No 



V/c : Zfc > no 



ifc = n - Tjino 
Vfc : ifc > 



E + E 



{v{2-^\...,2~^-)~vo) 



and consequently 



woiV| < + c (TV - TVo) < + c-^N = eiV. 



Division by N leads to (f^ . 

Asymptotically one gets instead of (|46|) for uq 7^ 



m — 1 

E] ^n+d-i-l 
i=0 



m — 1 



m — 1 



(-2r 



TO — 1 



i-oA^E 

n + TO — 1 
TO — 1 



vo 



(l-2P)"-\ 



in other words Ve > 3N Vn > iV 

m— 1 



E 

i=0 



^n+d— i— 1 



TO — 1 

i 



i~2Y 



n + TO — 1 
m — 1 



(2P- l)'"-i(l + e). 



The rest follows as in Theorem 



□ 
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